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Abstract 

The densest local packings of N identical nonoverlapping spheres within a radius R m i n {N) 
of a fixed central sphere of the same size are obtained using a nonlinear programming method 
operating in conjunction with a stochastic search of configuration space. Knowledge of R m i n (N) in 
d-dimensional Euclidean space R d allows for the construction both of a realizability condition for 
pair correlation functions of sphere packings and an upper bound on the maximal density of infinite 
sphere packings in R d . In this paper, we focus on the two-dimensional circular disk problem. We 
find and present the putative densest packings and corresponding R m i n (N) for selected values of N 
up to N = 348 and use this knowledge to construct such a realizability condition and upper bound. 
We additionally analyze the properties and characteristics of the maximally dense packings, finding 
significant variability in their symmetries and contact networks, and that the vast majority differ 
substantially from the triangular lattice even for large N. Our work has implications for packaging 
problems, nucleation theory, and surface physics. 

PACS numbers: 



I. INTRODUCTION 



A packing is denned as a set of nonoverlapping objects arranged in a space of given di- 
mension. Packings of identical nonoverlapping spheres in <i-dimensional Euclidean space M. d 
have been employed in condensed matter and materials physics as models for the structures 
of a diverse range of substances from crystals and colloids to liquids, amorphous solids and 
glasses [1-3]. In structural biology, molecular dynamics simulations of interactions between 
large numbers of molecules employ chains of nonoverlapping spheres as models for various 
biological structures such as proteins and lipids [4-6] . 

In part due to the ability of these conceptually simple models to describe many of the 
fundamental characteristics of more complex substances, understanding the properties of 
sphere packings has also long been an area of interest in mathematics (for example, see 
[7]). However, solving even some of the most basic of mathematical problems has proved 
challenging. For example, a proof of the Kepler conjecture, a proposition stating that the 
face-centered cubic lattice is the densest possible arrangement of spheres for d — 3, has 
only recently emerged [8]. Furthermore, the kissing number K d , or number of identical 
d-dimensional nonoverlapping spheres that can simultaneously be in contact with (kiss) a 
central sphere, was until recently only known rigorously for d = 1 — 3, 8 and 24 [9], though 
Musin [10] has now proved the d — 4 case (K4 = 24). 

One sphere packing problem that has not been generally addressed for an arbitrary num- 
ber of spheres is that of finding the maximally dense (optimal) packing(s) of N identical 
d-dimensional nonoverlapping spheres near (local to) an additional fixed central sphere such 
that the greatest radius R from any of the surrounding spheres' centers to the center of the 
fixed sphere is minimized. This problem is called the densest local packing (DLP) problem 
[11]. There is a single minimized greatest radius, denoted by R m i n (N), for each N in the 
DLP problem in R d , though generally for each N there may be multiple distinct packings 
that achieve this radius. Figure 1 depicts a conjectured optimal packing, belonging to point 
group D 5h [12], for the DLP problem for iV = 15, d = 2, with R min (15) = 1.873123 . . . [13]. 

In various limits, the densest local packing problem encompasses both the kissing number 
and (infinite) sphere packing problems. The former is a special case of the DLP problem in 
that Kd is equal to the greatest N for which R min (N) = 1, and the latter is equivalent to 
the DLP problem in the limit that N — > 00. The equivalence of the latter problem may be 
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FIG. 1: A conjectured DLP optimal packing (point group D^) for N = 15, d = 2, i? m j n (15) = 
1.873123 . . . , with encompassing sphere of radius i? m j„(15) + 0.5 = 2.373123 .... 

explained by observing that in the limit as N — > oo, the boundary of radius R m i n (N) — > oo, 
and that in this limit the ratio of the number of spheres within a fixed finite distance of the 
boundary to the number in the bulk is zero. 

The densest local packing problem is relevant to the realizability of functions that are 
candidates to be the pair correlation function of a packing of identical spheres. For a 
statistically homogeneous and isotropic packing, the pair correlation function is denoted 
<72(f)j it is proportional to the probability density of finding a separation r between any 
two sphere centers and normalized such that it takes the value of unity when no spatial 
correlations between centers are present. Specifically, no function can be the pair correlation 
function of a point process (where a packing of spheres of unit diameter is a point processes in 
which the minimum pair separation distance is unity) unless it meets certain necessary, but 
generally not sufficient, conditions known as realizability conditions [14-16]. Two of these 
conditions that appear to be particularly strong for the realizability of sphere packings [17] 
are the nonnegativity of g2(r) and its corresponding structure factor S(k), where 

S(k) = l+ph(k) (1) 
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with number density p and 

= [r^h(r) 1 0Q r dr (2) 

the (i-dimensional Fourier transform of the total correlation function h(k) = g 2 (r) — 1, with 
J v {x) the Bessel function of the first kind of order v. 

The (^-invariant process of Torquato and Stillinger [15] is a method to maximize the 
number density p associated with the structure factor S(k) of a given parameterized family 
of test functions, where each function in the family is a candidate to be the pair correlation 
function of a statistically homogeneous and isotropic packing of spheres. In the ^-invariant 
process, the problem of finding the maximal achievable p is posed as an optimization prob- 
lem: maximize p over the parameters subject to the nonnegativity of the test function and 
its corresponding structure factor. This process could be improved by the addition of further 
realizability conditions on the pair correlation function, assuming that these further condi- 
tions included information beyond that incorporated in the two nonnegativity conditions 
discussed above. 

Knowledge of the maximal number of sphere centers that may fit within radius R from 
an additional fixed sphere center, where that maximal number is equal to the greatest iV 
in the DLP problem for which R min (N) < R, may be employed to construct an additional 
realizability condition on g 2 (r). As was discussed in previous papers [11, 18], this realizability 
condition has been shown to encode information not included in the nonnegativity conditions 
on pair correlation functions and their corresponding structure factors alone. 

The DLP problem may be alternatively stated as the problem of finding the densest 
packing of iV identical nonoverlapping spheres of unit diameter near an additional fixed 
sphere, where number density p is measured over the volume enclosed by an encompassing 
sphere of radius R+0.5 (see Fig. 1) centered on the fixed sphere. We note that for N spheres 
of unit diameter, the number density p of the packing is linearly proportional (by a constant 
that varies only with dimension) to the packing fraction 0(i?+O.5), the fraction of the volume 
of the encompassing sphere covered by the spheres of unit diameter. As will be discussed 
in detail later, the maximal infinite-volume packing fraction of identical nonoverlapping 
spheres in d dimensions may be bounded from above by employing a specific definition of 
local packing fraction for a given number N of spheres. 
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For small numbers of spheres (N < 1000) in low dimensions (d < 10), a algorithm 
combining a nonlinear programming method with a stochastic search of configuration space 
can be employed on a personal computer to find solutions to the DLP problem. Using such 
an algorithm, the details of which are outlined in Appendix A, we find and present putative 
DLP optimal packings and their corresponding R m i n (N) in IR 2 for N = 1 to N = 109, and 
for the values of N corresponding to full shells of the triangular lattice from N = 120 to 
N = 348. Though we recognize that the putative optimal packings found by our algorithm 
are not rigorously proved to be optimal, we analyze each configuration of N spheres under 
the assumption that it is a global minimum of the DLP problem. This assumption of 
optimality is supported by the proved robustness of the algorithm in recovering the known 
and strongly conjectured global minima of the DLP problem (e.g., the kissing numbers for 
d — 1 — 4 and for d — 8, the curved hexagonal packings for d = 2; N = 18, 36, 60, 90 and 
126 [19]) and by repeated testing. 

The aforementioned realizability condition on the pair correlation function is valid 
whether or not the putative optimal packings we have found are indeed global minima, 
as global minima simply provide the most restrictive realizability condition. However, the 
upper bound on the maximal density of an infinite sphere packing requires knowledge of 
proved optimal R m i n (N) to be rigorously correct, though we have found in practice that 
for d — 2 and over the range of N tested that our putative bound is valid. With regard to 
this finding and the proved robustness of the algorithm over the range of N studied, in the 
following sections we refer to all DLP packings and R min (N) presented as optimal. 

In Sec. II, we discuss the realizability condition that results from knowledge of a finite 
number of R m i n (N) in a space M d of arbitrary dimension d. We present the condition derived 
from knowledge of R min (N) for JV = 1 to N = 109 in R 2 and compare the R m i n (N) values 
to the shell distances in a triangular lattice of N disks. In Sec. Ill, we construct a logical 
argument to prove the validity of the aforementioned upper bound, and we present the 
d — 2 upper bounds derived from our method for selected N from N = 6 to iV = 348. In 
Sec IV, we present d = 2 optimal packings and their corresponding R min (N) for selected 
values of N from N = 10 to N = 348. We analyze the optimal packings presented and 
discuss their symmetry characteristics, noting that there is significant variability in both 
the configurations and symmetry elements of optimal packings over the range of N studied. 
In Sec. V, we summarize our results and findings, and we discuss some of the implications 
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of our work. 

In a sequel to this paper, we will present and analyze DLP optimal packings and their 
corresponding R min (N) for the d = 3 case over a similar range of N. We will compare 
Rmin{N) values to shell distances in Barlow packings [20], where the Barlow packings, of 
which the best-known are face centered cubic (FCC) and hexagonal close packed (HCP) 
arrangements, all individually achieve the maximal infinite-volume packing fraction for d = 
3, cf)f = n/Vl8 = 0.740481 .... 

The coordinates for and images of DLP optimal packings and values for R rnin (N) over 
the entire range of N studied can be found on the authors' website [21]. 



II. PAIR CORRELATION FUNCTION REALIZ ABILITY AND Z max (R) 

The realizability condition on (fetj) results from a relation between an upper bound on 
the maximal value of the function Z(R), to be defined shortly, and g2{r). The function 
Z(ri, R) is defined for packings of nonoverlapping spheres of unit diameter as the number 
of sphere centers that are within distance R from a (additional) sphere center at position 
r-j, with % an index over centers. The maximum over all rj of Z{v^ R) is an upper bound 
on the maximum of the (different) function Z(R), where Z(R) is defined for a statistically 
homogeneous packing as the expected number of sphere centers within distance R from any 
given sphere center, or equivalently as the average of Z(rj, R) over all i. The function Z(R) 
can be related to the pair correlation function g2(r), where for a pair correlation function 
g 2 (r) that is direction-dependent, g 2 (r) is the directional average of #2(1"), by 



r R 

Z(R) = psi(l) / x d - 1 g 2 {x)dx. 
Jo 



(3) 



In Eq. (3), p is the constant number density of sphere centers and si(r) is the surface area 
of a sphere of radius r in M. d , 



27T d/2 r 



d-1 



* l(r) = -ww (4) 

The maximum at fixed R of the function Z(r{, R) over all possible configurations of sphere 
centers {rj} is equal to the greatest number N for which R m i n (N) < R for a DLP optimal 
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packing of N spheres in M. d . Defining Z max (R) for all R as this greatest N, it follows that 

Z{R) < Z max {R) (5) 

for any sphere packing. 

Equation (5) is a realizability condition on g<i (r), with Z(R) in IR d defined in terms of 
g2{r) in Eq. (3) and Z max (R) defined completely by the solutions to the DLP problem over 
all N. In IR 2 , the function Z max (R) may be compared to the function Z tri (R), with Z tri (R) 
defined as the sum of the number of disk centers included in all (full) shells of radius less 
than or equal to R in a triangular lattice of contacting disks. Both Z max (R) and Z tr i{R) 
increase roughly linearly with R 2 , as the area of a disk of radius R is proportional to R 2 . 
Clearly, Z max {R) > Z tr i{R) for all R, as can be seen in Fig. 2, a plot of Z max {R) vs. R 2 for 
N = 1 — 109 and d = 2 alongside a plot of Z tri (R). 




FIG. 2: Z rnax (R) vs i? 2 , as determined by optimal and putative optimal solutions to the DLP 
problem for N = 1 to N = 109, and Zt r i{R). The radius R of the (larger) disk enclosing the 
centers of the N (smaller) disks and fixed disk is measured in units of the diameter of the enclosed 
disks. 

In arbitrary dimension d, the function Z max (R) is zero for R < 1 due to the nonoverlap 
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condition. For R — 1, Z max (R) in IR d is equal to the kissing number K d . For R > 1, 
Zmax(R) should grow approximately as _R d in proportion with the growth of the volume of 
a d-dimensional sphere, though in a separate work [11] we have proved that this cannot be 



dimension cannot exceed the maximal number of sphere centers that can be placed on the 
surface of a sphere of radius R. Alternatively stated, this counterintuitive result requires 
that for R < r, Z max (R) can grow only as the surface area R d ~ l . Specifically for d — 2, 3 
and 4, Z max (R < r) is less than or equal to 10, 33 and 120, respectively. 

III. BOUNDS ON INFINITE SPHERE PACKINGS AND THE DLP PROBLEM 

We discuss two distinct methods through which the function Z max {R) in R d can be 
employed to bound from above the maximal infinite- volume packing fraction of an infinite 
packing of identical nonoverlapping spheres. The first has been discussed in detail in two 
separate works [11, 18]; it is precisely the method of Cohn and Elkies in [22]. In Ref. [22], 
the authors employ an infinite-dimensional linear program that is the dual of the ^-invariant 
program [17] discussed in Sec. I to find the best known bounds on the maximal infinite- 
volume packing fraction of sphere packings at least in dimensions four through 36. An 
improved method to bound from above adds the information encoded in the Z max {R) 
realizability condition to augment the approach of Cohn and Elkies in [22] as proposed by 
Cohn, Kumar and Torquato [23]. 

The second method bounds from above by the maximal local packing fraction 0*(A) 
of a packing of a number A of identical nonoverlapping spheres around an additional central 
sphere. The local packing fraction 0(A), of which 0*(A) is the maximum, is defined for A 
spheres around an additional fixed central sphere as the total volume of the A + 1 spheres 
divided by the volume of a sphere of radius R, where R is, as in the DLP problem, the 
greatest of the distances from the centers of the A surrounding spheres to the center of the 
fixed sphere. From this definition of 0(A), the maximal local packing fraction 0*(A) for A 
d-dimensional spheres of unit diameter takes the form 



the case for R < r, with r 



(1 + a/I 5 ))/ 2 the golden ratio. For R < r, Z, 



max 



(R) in any 



</>*(A) 



A+ 1 



(6) 



(2R mm (N)) d ' 
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TABLE I: Packing fraction terms and definitions used in the text. 



Symbol 


Term 


Definition 


<j>(R) 


packing fraction 


volume fraction of a (larger) sphere of radius R covered 


by identical nonoverlapping spheres of unit diameter 


(Pmnr -K 

Y • > vxiiJLi \ J 


maximal 
packing fraction 


greatest achievable 4>{R) for a given R 


0°° 


infinite- volume 


fraction of space covered by identical nonoverlapping 


packing fraction 


spheres in a given infinite packing 




maximal 




#° 


infinite- volume 
packing fraction 


greatest achievable infinite-volume packing fraction 






ratio of the sum of the volumes of a nonoverlapping fixed 


4>{N) 


local 


central sphere of unit diameter and its surrounding N 


packing fraction 


same-size spheres to the volume of a (larger) sphere of 
radius R, with R defined as in the DLP problem 




maximal local 
packing fraction 


greatest achievable 4>(N) for a given N in R d ; see Eq. (6) 




average local 


average packing fraction within a window of radius R of 


packing fraction 


identical nonoverlapping spheres of unit diameter 



where R m i n (N) is as before the optimal radius in the DLP problem for N spheres in M d . For 
the sake of convenience, we have collected and defined the various packing fraction terms 
used in this section in table I. 

The statement that 0* ( N) bounds from above for certain N relies on a construction 
that links local packing fraction 4>{N) to the infinite-volume packing fraction 0°° of a packing 
of identical nonoverlapping spheres. The construction proceeds as follows. First, a spherical 
window of radius R min (N) is centered on an arbitrary sphere in a single configuration of an 
infinite packing of identical nonoverlapping spheres of unit diameter. An infinite packing 
of identical nonoverlapping spherical windows is created by replicating the initial window 
infinitely many times and placing the (replicated) window centers in the exact (scaled) spatial 
configuration of the centers of the original infinite packing of spheres of unit diameter. The 
only difference between the two configurations is that the configuration of windows is scaled 
by 2R min (N), the ratio of the radius of a window to the radius of a sphere of unit diameter. 

As will be made precise in the following paragraphs, for any such packing of windows 
and spheres of unit diameter, a rigid rotation for the overlayed packing of windows can be 
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found such that the average local packing fraction <f>(R min (N)) of spheres of unit diameter 
within the windows is equal to the infinite-volume packing fraction 0°° of the spheres of 
unit diameter in R d . The concepts of overlay and rotation are illustrated in Fig. 3 for a 
triangular lattice of disks of unit diameter. 




FIG. 3: Illustration of a rotation of an infinite packing of identical nonoverlapping windows of 
radius R, arranged on the sites of a triangular lattice, overlayed upon an infinite packing of smaller 
identical nonoverlapping disks. A rotation is selected such that irrational ratios are achieved 
between the components of at least one of the lattice vectors of the packing of windows in the 
directions of the lattice vectors of the packing of smaller disks. As a result, at large distances from 
the axis of rotation, any window can be thought of as being placed at random onto the packing of 
smaller disks. It follows that the average local packing fraction 4>(R) of the smaller disks within 
the windows is equal to the infinite- volume packing fraction (f)°° of the smaller disks. 

It suffices to apply the aforementioned construction to periodic packings, as it has been 
shown that periodic packings in M. d can obtain an infinite-volume packing fraction 0°° ar- 
bitrarily close to the maximal infinite- volume packing fraction (for example, see [22]). 
A periodic packing can be defined in terms of a lattice A, where A in R d is a subgroup 
consisting of the integer linear combinations of a set of vectors that constitute a basis for 
M. d . For identical nonoverlapping spheres, a lattice packing is a packing where the centers 
of the spheres are located at the points of A. In such a lattice packing, the space R d can 
be divided into finite-size identical nonoverlapping regions called fundamental cells, each 
containing the center of only one sphere. 

A periodic packing is a more general formulation of a lattice packing. For identical 
nonoverlapping spheres, a periodic packing is obtained by placing a fixed configuration of 
a number M of spheres in a fundamental cell that is then periodically replicated (without 
overlap between cells or spheres) to cover R d . The fixed configuration of M spheres within 
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each fundamental cell is arbitrary subject only to the overall nonoverlap condition of the 
periodic packing of spheres. As used here, the term "lattice" is the same as "Bravais lattice" 
conventionally used in the physics literature. 

Consider an infinite periodic packing (with lattice basis vectors {uj}) of nonoverlapping 
spheres of unit diameter in IR d , d > 1, with infinite- volume packing fraction 0°°. Place 
an infinite periodic packing (with lattice basis vectors {v.,}) of identical nonoverlapping 
windows of radius R min (N) over the infinite periodic packing of spheres of unit diameter in 
the manner of the construction discussed previously. For the radius R min (N) of the windows, 
any positive integer iV* G N can be considered, where in R d the set N is defined such that 
each iV* is the greatest number N of spheres for which R m i n {N*) = R m i n (N). For example, 
with iV = 1 in two dimensions, R m i n (l) = 1, and the greatest N for which R m i n (N) = 1 is 
N = iV* = 6 [24]. It is intuitively clear that for any greatest such N = N* and R = R m in(N*) 
that 

0*(iV*) > <pmax(Rmin(N*)), (7) 

where (p ma x(R) is defined as the maximal fraction of space that identical nonoverlapping 
spheres of unit diameter may cover in a spherical window of radius R [25]. 

Now rigidly rotate the packing of windows around the center of the original window as 
in Fig. 3. The component of a window-packing basis vector v n that is in the direction of 
a given unit-diameter sphere packing basis vector u m is v n ■ u m /|u m |. A rotation may be 
found such that the ratios {v n • Uj/|u,'| 2 } of the components of (at least) one of the window 
basis vectors v„ in the directions of each of the sphere packing basis vectors {uj} to the 
respective magnitudes {|u.,|} of the sphere packing basis vectors are all irrational. This 
concept is illustrated in Fig. 4 in M 2 for parallelogram fundamental cells. 

After the rotation, due to the irrationality of the ratios of lattice vector components, 
the average fraction of space (f>(R m i n (N)) covered by the spheres of unit diameter in each 
window of the window packing is equal to the infinite-volume packing fraction 0°° of the 
unit-diameter spheres. 

With this construction, the packing fraction 0°° of the spheres of unit diameter, equal to 
the average local packing fraction (f)(R min (N)) of the unit-diameter spheres within a window, 
can be compared to the maximal local packing fraction 0*(iV) of that same window. As 
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FIG. 4: Illustration of a rigid rotation of a window packing with parallelogram fundamental cells 
such that the ratios of the components of vi in the directions of each of ui and 112 to the magnitudes 
of ui and 112, respectively, are both irrational. 
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FIG. 5: maximal packing fraction = 0.906900 . . . compared to putative upper bounds on 
for d = 2 as determined by putative optimal solutions to the DLP problem for selected JV» £ N up 
to iV* = 348. The minimum upper bound determined here is at = 336, 0*(336) = 0.928114. 

(pmax(Rmin(N)) is clearly greater than or equal to 0°°, using Eq. (7) we have 



M N *) > <Pmax{Rmin{N*)) > <p C 



(8) 



for JV„ E N in R d . 

Equation (8) is true for all feasible 0°° including the maximal infinite-volume packing 
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fraction 0^° of nonoverlapping spheres in R d . The maximal local packing fraction 0*(-/V*) is 
therefore an upper bound on the maximal infinite- volume packing fraction 0£°, or 

0*(iV*)>C N*eN. (9) 

Figure 5 plots for d = 2 and for selected values of iV* e N up to = 348 the putative upper 
bound defined by ^(iV*) alongside the proved maximal packing fraction, 0^° = 7r/vT2. 

While the maximal infinite- volume packing fractions for (/-dimensional identical nonover- 
lapping spheres are known with analytical rigor for d = 1 — 3, d — 8 and d = 24 [26], this 
method could be used to improve upon the upper bounds on maximal packing fraction 
in dimensions where a value for has not not yet been proved. It is important to reiterate 
though that to be rigorous, the upper bound as defined requires a DLP radius R min (N) that 
is proved optimal. 

IV. OPTIMAL PACKINGS IN TWO DIMENSIONS 

The packing of nonoverlapping disks that uniquely achieves the highest infinite-volume 
packing fraction 0£° = 7t/a/T2 for d = 2 is well-known to be the packing such that each disk 
is in contact with exactly six others, with centers arranged on the sites of the triangular 
lattice. As DLP optimal packings are packing fraction-maximizing arrangements of N disks 
centered on a (additional) disk, one might expect that all or a major subset of disks in any 
given DLP optimal packing will always sit on the sites of the triangular lattice. 

Over the range of N studied, however, this is in fact very infrequently the case. We 
find that at only three values of N greater than 6, at N — 12, 30 and 54, are triangular 
lattice configurations also DLP optimal configurations, while most DLP optimal packings 
are significantly more locally dense than a packing of N disks around a central disk with 
centers on the sites of the triangular lattice, as is illustrated in Fig. 2 for N = 1 — 109. 

In general, we find wide variation in the symmetries and other characteristics of DLP 
optimal packings. For the majority of cases, there appear to be an uncountably infinite 
number (a continuum) of optimal configurations of N sphere centers at optimal R m i n (N), 
with the continuum attributable to the presence of rattlers. A rattler in a packing of spheres 
in M. d is a sphere that is positioned such that it may be individually moved in at least one 
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direction without resulting overlap of any other sphere within the packing or the packing 
boundary (in this case, the encompassing sphere of radius R min (N) + 0.5), i.e., a rattler is a 
sphere that is not locally jammed [27, 28]. The rattlers present in the following figures are 
indicated by a lighter shading, while the fixed central spheres (disks) are indicated by an 
open circle. 

A. Packings that are proved optimal 

In two dimensions for iV < 10, we proved [11] that optimal packings are those in which 
Rmin(N) is equal to the radius of the minimal-radius circle onto which the centers of N disks 
may be packed. These R min (N) may be analytically calculated via simple trigonometry, 
yielding R min (N) = 1 for iV < 6, and 

R min (N) = 1 6 < N < 10. (10) 

y/2(l — cos27i/N) 

For 6 < N < 9, DLP optimal packings with R min (N) as defined in Eq. (10) are unique 
up to rotations and correspond to configurations where all sphere centers lie on the circle of 
radius R m i n (N) at distance unity from two adjacent centers. In general, for many N there 
is a unique (up to rotations) optimal packing, and in a smaller number of cases such as for 
N = 10, 60, 90 and 126, there are a finite number of degenerate optimal packings. Figure 6 
depicts three of a finite number of optimal packings for N — 10, d — 2, where i? m j n (10) = r; 
these packings are formed by radially translating up to five disks with centers on the circle 
at radius i? m j„(10) = r, the golden ratio, to be in contact with the fixed central disk at 
distance unity. 

B. Curved hexagonal packings 

For N > 10, we know of no packings that have been proved optimal. However, for certain 
N, previous studies have found packings that we conjecture to be DLP optimal packings. For 
N + 1 equal to hexagonal number 3k(k + 1) + 1, k > 1 an integer, Lubachevsky and Graham 
[19] found a class of packings called curved hexagonal packings that they conjectured to be 
the densest packings of N + 1 identical nonoverlapping disks within an encompassing disk 
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(c) 

FIG. 6: Three of a finite number of optimal cases for N = 10, i? m j n (10) = r = 1.618204 . . . formed 
by radially translating disks whose centers lie on a circle of radius R m i n (W) = r to be in contact 
with the fixed central disk, (a) point group C^- (b) point group Civ- (c) point group D^. 

for k = 1 ... 5. The characteristics of this class include that each packing has a disk fixed 
at the center of the encompassing disk, and that all curved hexagonal packings belong to 
point group C^, meaning that they are invariant under 60° rotation or inversion through 
the origin. Further, for k > 4, for each k there are a finite number of degenerate packings of 
equal density (for 1 < k < 3, there is a unique densest packing). The degenerate packings 
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are chiral; each ring of disks, or disks that share a common radius, beginning with the 
fourth ring (the fourth farthest from the center), can be angularly oriented in more than one 
distinct fashion relative to the preceding ring. By reorienting rings, the degenerate packings 
for a given k can be generated from one another. Figure 7 depicts curved hexagonal packings 
for N = 60, 90 and 126. 

DLP optimal packings for N disks are equivalent to the densest packings of N + 1 disks 
enclosed in an encompassing disk when one of the disks is fixed at the center. For N = Q, 
18, 36, 60 and 90 (k — 1 ... 5), we find that the curved hexagonal packings are the only DLP 
optimal packings, in support of the conjecture of Lubachevsky and Graham that curved 
hexagonal packings are the densest packings up to k = 5. Further, for N = 126 (k — 6), we 
also find that curved hexagonal packings are the only DLP optimal packings, indicating that 
though there are packings of 127 unconstrained disks within an encompassing disk denser 
than the curved hexagonal packings (as were found by Lubachevsky and Graham [19]), the 
curved hexagonal packings remain the densest packings of 126 disks around a fixed central 
disk. This is not the case for N = 168 (k — 7), as we find DLP optimal packings with higher 
density, such as the N = 168 packing to be shown later in the top panel of Fig. 16. 

C. Wedge hexagonal packings 

Another class of packings, previously unidentified, contain a subset of disks with centers 
arranged on the sites of the triangular lattice and the remainder arranged in six "wedges" . 
We hereafter term such packings wedge hexagonal packings. Wedge hexagonal packings are 
not DLP optimal packings when arranged symmetrically (point group D^); however, minor 
deviations from perfect symmetry in a wedge hexagonal packing can produce a DLP optimal 
packing. Figure 8 depicts such DLP optimal packings for N = 84, 120 and 162. Lines to 
guide the eye have been drawn on the three optimal packings in Fig. 8. 

In a wedge hexagonal packing, the subset of disks with centers arranged on the sites 
of the triangular lattice contains two parts; a regular hexagonal core of hexagonal number 
3k(k + 1) + 1 disks, with k > 3 odd; and six 'branches' composed of (pk — a) disks, with 
p > 2 and a > 1 integers, extending from each of the vertices of the core regular hexagon. 
The branches are k disks wide and p disks long, with a of the farthest disks removed such 
that the end of the branch approximates a circle (as opposed to the point of a triangle). The 
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(a) 




(c) 



FIG. 7: Curved hexagonal packings that are also DLP optimal packings, (a) N = 60, R m i n (60) = 
3.830649..., point group C 6h . (b) N = 90, R min (90) = 4.783386..., point group C 6h . (c) 
N = 126, R min {126) = 5.736857. . . , point group C 6h . 

six wedges are arranged roughly as p(p+ l)/2 bowling pins and lie in between the branches, 
with each of the six 'lead pins' placed at the midpoint of each side of the core hexagon. The 
DLP optimal packings in Fig. 8 are the wedge hexagonal packings, with minor deviations 
in the positions of some disks, for (k,p, a) = (3, 2 — 4, 1). 



17 



(a) 




(c) 

FIG. 8: DLP optimal packings that are minor deviations from "wedge hexagonal" packings, (a) 
N = 84, i? mm (84) = 4.581556. . . , (k,p,a) = (3,2, 1), point group G t . (b) N = 120, ^(120) = 
5.562401... , (k,p,a) = (3,3,1), point group Q. (c) N = 162, R min (162) = 6.539939... , (k,p,a) = 
(3,4,1), point group Cj. 

In general, the deviations necessary to produce a DLP optimal packings from a wedge 
hexagonal packing occur in the branches, and to a lesser degree, the wedges of the packing, 
while the core regular hexagon retains perfect six-fold symmetry. The deviations required 
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differ for each wedge hexagonal packing, but from our observations produce packings where 
the backbone maintains inversion symmetry through the origin, where the backbone of a 
packing is defined as the packing excluding the rattlers. Such deviations can be seen in the 
branches and wedges in the DLP optimal packings for N = 198 and N = 312 in Fig. 9, 
which correspond to the (slightly altered) wedge hexagonal packings with (k,p, a) = (5, 3, 3) 
and (7,3,3), respectively. 




(b) 

FIG. 9: DLP optimal packings that are minor deviations from the "wedge hexagonal" packings, 
(a) N = 198, d = 2. R min (198) = 7.201130 (k,p,a) = (5, 3, 3), point group d. (b) N = 312, 
Rmin(312) = 9.141107..., (k,p,a) = (7,3,3), point group C t . 

D. DLP optimal packings with high symmetry 

Many DLP optimal packings exhibit symmetries other than inversion symmetry through 
the origin as exhibited by the altered wedge hexagonal packings shown in Figs. 8 and 
9. These symmetries include perfect bond orientational order, invariance under rotation 
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through an angle, and invariance under reflection across an axis. A list of packing point 
group, alongside other packing properties such as R min (N) value, of all of the DLP optimal 
packings depicted in this work appears in Appendix B. 




(c) 



FIG. 10: Three examples of DLP optimal packings incorporating interesting symmetry elements, 
(a) TV = 25, R mm {2h) = 2.497212..., point group D 5h . (b) N = 11, R min {ll) = 1.685854..., 
point group C 2v . (c) N = 32, i? min (32) = 2.794164 . . . , point group D 2h . 
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Perfect five-fold symmetry, disallowed to regular infinite crystals, is exhibited by three of 
the optimal packings studied. Five-fold rotational symmetry is present in the N — 15 (Fig. 
1), TV = 10 (bottom panel of Fig. 6) and N = 25 (top panel of Fig. 10) packings. The 
N = 25 optimal packing also has perfect five-fold bond orientational order, evident in that all 
nearest-neighbor disk pairs are at one of five angles relative to any fixed coordinate system. 
Additionally, it is of note that the N = 25 packing may be tiled by 15 identical rhombuses 
of acute angle 72° with vertices placed at disk centers, where the rhombus of acute angle 
72° is known to be the 'thicker' of the two types of rhombus present in a Penrose tiling [29] . 
The top panel of Fig. 14, a diagram of the contact network for the N = 25 optimal packing, 
depicts these rhombuses. 

Two other optimal packings incorporating interesting symmetry elements are the N — 1 1 
and N = 32 packings (center and bottom panels of Fig. 10, respectively). The packing 
depicted for iV = 11 belongs to symmetry group C 2v and appears to be the unique DLP 
optimal packing of 11 disks. The packing for N = 32 is one of an infinity of possible 
packings due to the presence of two rattlers; however, the backbone of the packing has 
reflection symmetry across two axes and inversion symmetry through the origin and belongs 
to symmetry group D 2 h- 

Though five-fold symmetry may be limited to packings with small N, high symmetry 
in general is not. For example, the backbone of the optimal packing with the largest iV 
presented here, N = 348, has six-fold rotation symmetry and belongs to point group D 6h . 
Figure 11 depicts the N = 348 optimal packing, which contains 24 rattlers. 



FIG. 11: DLP optimal packing for N = 348, R min (M8) = 9.620709. . . , belonging to point group 
D 6h . 
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E. Unusual features in select optimal packings 




(c) 

FIG. 12: Three examples of DLP optimal packings with circular cavities about the fixed central 
disk, (a) TV = 24, i? min (24) = 2.425256. . . , point group D 3h . (b) N = 45, -R mi „(45) = 3.374023, 
point group C\. (c) N = 95, i? m j n (95) = 4.958096, point group C\. 

A prevalent feature found in many of the DLP optimal packings studied is a cavity 
consisting of a ring of disks enclosing, but not contacting, the fixed central disk. The 
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counterintuitive presence of this feature is related to the aforementioned fact that in any 
dimension d, no spherical window of radius R < r centered on a central nonoverlapping 
sphere of unit diameter may encircle more sphere centers than the number (plus one) that 
can be placed on the encircling sphere's surface. Circular cavities around the fixed central 
disk appear in many DLP optimal packings, including those for N = 24, 45 and 95 disks, 
as illustrated in Fig. 12. There are 9, 8 and 7 disks, respectively, forming the walls of the 
cavities in the three packings in Fig. 12. 

Two particularly notable DLP optimal packings that include a cavity enclosing the fixed 
central disk are the N = 40 and iV = 66 packings shown in Fig. 13. These packings are 




(b) 

FIG. 13: DLP optimal packings consisting of layers of distorted rings, or 'eyes', enclosing the central 
disk, (a) N = 40, i? min (40) = 3.136712 . . . , point group D 2h - (b) N = 66, R min {66) = 4.104997, 
point group C\. 

composed of layers of distorted rings, where the distorted rings appear as eye-like closed 
curves of varying curvature with each successive layer from the center more circular than 
the last. It is curious that even though the only shape imposed upon the packings, in the 
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form of the encompassing disk and the disks themselves, is circular, that an optimal packing 
incorporating distorted rings emerges for these numbers of disks. The bottom panel of Fig. 
14 is a diagram of the contact network for the backbone of the N = 40 optimal packing. 




(b) 

FIG. 14: Diagrams [30] of contact networks for the N = 25 (a) and N = 40 (b) optimal packings, 
with point group symmetries D$h and -D2/1, respectively. 

The presence of these cavities about the central disk leads to an interesting counterintu- 
itive result. Suppose that in a binary liquid of nonoverlapping disks of unit diameter, one 
species of disk is endowed with an attractive square well potential. The potential of the 
"attractive" disks, present in the dilute limit in comparison to the "nonattractive" disks, 
acts on the centers of the "nonattractive" disks extending only to a distance just larger than 
any of the R min (N) presented in Figs. 12 or 13. If the depth (strength) of this square well 
were made arbitrarily large, then the result would seem paradoxical: unbounded attraction 
to the disks, in a minimal energy configuration, would eliminate contact with these disks. 

This effect in such a binary liquid of disks requires that the pair correlation function 
depicting the probability density of finding the centers of a given number of "nonattractive" 
disks a certain distance from the centers of "attractive" disks be zero from r = to r = 
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Ro > 1 a specified distance in excess of the diameter of the disks. None of the currently 
available pair correlation function theories are able to predict this effect, including the crucial 
dependence of R on R min (N), because the underlying approximations of the currently 
available theories cannot account for the basic many-body geometrical features involved. 

F. Imperfect symmetry 

Not all DLP optimal packings exhibit perfect symmetry; for many N, a subset of disks in 
an optimal configuration appear to mimic a symmetric packing, but the packing as a whole 
exhibits only imperfect symmetry. One situation in which this occurs frequently is when 
the number of disks N in the optimal packing is close to a different number for which the 
optimal packing is relatively unusually dense. For example, the optimal N = 59 packing 
shown in the top panel of of Fig. 15 lacks any of the symmetry elements described above 
but nonetheless closely resembles the particularly dense N = 60 curved hexagonal packing 
(top panel of Fig. 7). 

Other packings in which imperfect symmetry is present include the N = 80 packing 
shown in the center panel of Fig. 15, the disks closer to the center of which are ordered 
with centers on the sites of a slightly distorted triangular lattice, and the N = 46 packing 
shown in the bottom panel of Fig. 15, which has imperfect five-fold symmetry. The N = 46 
packing, along with the N = 45 packing (center panel of Fig. 12), together illustrate another 
finding: that the structure of optimal packings even for consecutive numbers of disks can 
vary substantially. 

G. Surface effects 

DLP optimal packings with N in the higher range of the packings studied appear, as 
N increases, to more and more resemble the triangular lattice in the bulk of the packing. 
Nonetheless, the surface of the packing always deviates significantly from the bulk crystal. 
In general, the optimal packings at higher N consist of a "bulk zone" with disk centers 
arranged in the triangular lattice surrounded by a "surface zone" with disk centers arranged 
in circular rings. This effect can be seen in all of the optimal packings with N in the higher 
range of N studied, including in those shown in Figs. 8, 9, 11, in the center panel of Fig. 
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(a) 




(c) 



FIG. 15: DLP optimal packings that exhibit imperfect symmetry, (a) N = 59, -R m m(59) = 
4.824374, point group C\. (b) N = 80, R min (80) = 4.514170, point group C\. (c) N = 46, 
R min (46) = 3.414304, point group C\. 

15 and in Fig. 16. 

Qualitatively, it appears that the radial width of the surface zone tends to increase with 
number of disks, though not as fast as the radial width of the bulk. Due to computational 
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time constraints, we were unfortunately not able to quantitatively verify this result at much 
larger N; however, should the observed trend continue, the width of the surface zone would 
continue to grow as the bulk grows, eventually becoming infinitely large as N — > oo. This 
does not imply that at very large N the surface zone would represent a substantial fraction 
of the total packing; as iV — > oo, the ratio of space covered by the surface zone to the space 
covered by the bulk zone will still be zero. 

For certain optimal packings exhibiting this surface effect, we also find that the optimal 
radius R min (N) is significantly smaller than the radius of the smallest disk centered on a 
fixed central disk that can enclose (an additional) N disks with centers arranged on the sites 
of the triangular lattice. Figure 16 depicts three optimal packings for which this difference 
is relatively large, at N = 168, 264 and 270, where the R min {N) differ by 0.248190, 0.128616 
and 0.162510 from the respective radii for the triangular lattice configuration, 7, v^73 and 
\/75. Each of the packings in Fig. 16 also displays its own interesting features, including 
a close resemblance on the left of the image of the N = 168 packing to the (3,4, 1) wedge 
hexagonal packing, and imperfect three-fold rotational symmetry for the N = 264 and 
N = 270 packings. 

V. CONCLUSIONS AND DISCUSSION 

The DLP problem is a local packing problem that in certain limits encompasses both the 
infinite sphere packing and kissing number problems. DLP optimal packings exhibit a wide 
variety of symmetries, vary significantly from packings with disk centers placed on the sites 
of a triangular lattice, and often vary significantly with consecutive N. Two local packing 
classes, the curved hexagonal and wedge hexagonal packings, lead to the densest or very 
dense packings of disks over the range of N studied. 

The optimal radii R m i n (N) corresponding to a packing of iV spheres in M. d form a realiz- 
ability condition on functions that are candidates to be the pair correlation function g2{r) 
of a statistically homogeneous and isotropic packing of spheres. This realizability condition 
incorporates more information than is included in the structure factor and pair correlation 
function nonnegativity conditions alone. Though the condition discussed here only applies 
to packings of identical spheres, equivalent Z max (R) functions for packings of differentiated 
non-spherical objects can be found and corresponding realizability conditions imposed. 
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(a) 




(c) 

FIG. 16: Three packings for which R m i n {N) is significantly less than the radius of a disk enclosing 
the centers of N + 1 disks arranged on the sites of the triangular lattice, (a) ./V = 168, i? m j n (168) = 
6.680013, point group C x . (b) N = 264, R min (26A) = 8.417769, point group C x . (c) N = 270, 
R min (270) = 8.497744, point group d. 

The function Z max (R) can also be employed in the two ways discussed to bound from 
above the packing fraction of an infinite sphere packing in dimension d. Similarly, Z max (R) 
for differentiated non-spherical objects can be employed to form upper bounds on corre- 



28 



sponding infinite packing maximal packing fractions. 

Our work has direct application to packaging, particularly to problems involving identical 
nonoverlapping disks within a circular boundary. Further, as is discussed in Appendix A, the 
algorithm we have employed to find the d = 2 putative DLP optimal packings presented in 
this work may be modified to study dense packings of ci-dimensional differentiated objects 
of various shape within different boundaries. Additionally, our work has implications for 
nucleation theory and surface physics, particularly in terms of the effects of imposing a 
circular boundary upon a packing of spheres with centers initially placed on the sites of a 
triangular lattice. In future work, we expect to investigate these implications and others in 
more depth. 

In a sequel to this paper, we will present and analyze DLP optimal packings and their 
corresponding R min (N) for three-dimensional spheres over a larger range of N. We will 
catalogue optimal packings of particularly high packing fraction and of unusual symmetry, 
and we will investigate the possibility of d — 3 extensions to special d = 2 classes of 
packings such as the curved hexagonal and wedge hexagonal packings. We will compare 
Rmin(N) values to shell distances in Barlow packings [20], where the Barlow packings, of 
which the best-known are face centered cubic (FCC) and hexagonal close packed (HCP) 
arrangements, all individually achieve the maximal infinite-volume packing fraction for d = 
3, (f>? = = 0.740481 .... 

Packing in IR 3 is intrinsically more complicated than in IR 2 . In IR 3 , there is a wider range 
of possibilities for contact coordination, particularly for twelve contacting spheres around 
a central sphere (with K 3 = 12) where there is an infinity of possible configurations, as 
compared to just one in IR 2 . Additionally, there is a single optimal infinite- volume packing 
configuration in IR 2 , i.e., the triangular lattice, whereas there is an infinite number in IR 3 , 
i.e., the Barlow packings. 

Preliminary findings indicate that there is less symmetry (as quantified by point groups) 
present for the majority of DLP optimal packings in IR 3 and more variation with N. Further, 
observations suggest that for the same value of N that there are significantly more locally 
jammed packing configurations in IR 3 with radius R > R m i n (N) than is the case in IR 2 . This 
finding has implications for the dynamics of nucleation occuring in pure supersaturated 
liquids. If a nucleus in a supersaturated liquid of identical nonoverlapping spheres in IR 3 is 
approximated as a group of N densely-packed spheres with centers within distance R of a 
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central sphere, then there are more packing configurations available to a nucleus of radius 
R, with R confined to a small finite range R > R min (N), in IR 3 than in IR 2 . 
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Appendix A: Description of the algorithm 

The algorithm employed to find the DLP optimal packings presented operates in two 
steps repeated iteratively in succession. In the first step, we use a method in nonlinear 
programming often called an augmented Lagrangian method [31] to find a local minimum of 
a specially formulated problem in Nd + 1 variables, where the first Nd variables correspond 
to coordinates of the centers of A spheres in d dimensions. In the second step, the configu- 
ration of A spheres found in the first step is spatially repositioned using a random number 
generator. Keeping track of the least local minimum, the two steps are repeated iteratively 
until (presumably) no further improvement on the least minimum can be achieved [32]. 

The nonlinear programming algorithm used in the first step seeks to solve a problem 
posed for Nd + 1 variables with 2A + jvC 2 inequality constraints, where ^C 2 is shorthand 
for the standard combinatorial formula. Calling the first Nd variables with % = 1 . . . N 
and k — 1 . . . d and the last variable u, the problem may be posed as follows: 

minw s.t. 

d 

k=l 

d 

-u + J2 x2 ik < V i : 1 ... A 

d 

1 - J2( Xik ~ x ^ ^ V 1 ... A, i < j. (Al) 

k=i 

The constraints in (Al) are split into three categories; the first category, containing A 
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constraints, requires that none of the centers of the spheres move within distance unity of 
the origin, i.e., that none of the N spheres overlap the fixed central sphere. The second 
category, also containing N constraints, requires that the difference between the squared 
distance from the origin to the centers of each of the N spheres and the objective function, 
simply the independent variable u, be less than or equal to zero. Since u is minimized, this 
condition sets oj equal to R 2 , i.e., the greatest of the squared distances from the origin to 
the any of the N sphere centers. The third category, containing ^C 2 constraints, requires 
that none of the N spheres overlap. 

An augmented Lagrangian (AL) method is employed to attempt solutions to the problem 
Al. In brief, an AL method is an iterative process designed to minimize a function that 
is the Lagrangian of (Al) augmented with a quadratic penalty function. The augmented 
Lagrangian is written, 

m 

L 7 (x, \) =U + ^Y1 
1=1 

where the {A;} are the Lagrange multipliers, 7 is the penalty parameter (often denoted by 
p), and the q(x) (often denoted by <?z(x)) are the m = 2N + ^C 2 inequality constraints. 
The function (A2) is minimized iteratively over x, yielding x p a local minimum of £ 7 (x, X p ), 
with \ p fixed. For each successive x p that is found, a new estimate for X p is made based 
on the violated constraints, i.e., the constraints in (Al) with positive values. Concurrently, 
the penalty parameter 7 is increased by a pre-specified multiple if the cumulative squared 
violation, or total error, is not a set amount smaller than the previous total error. Eventually, 
7 reaches a value such that the total error is smaller than a specified (small) tolerance, at 
which point the algorithm terminates. 

In the version of the AL method we use to find putative DLP optimal packings, the 
function L 7 (x, A) in each iteration is minimized using a conjugate gradient method with a 
directional minimizing algorithm that employs cubic interpolation. Though the conjugate 
gradient method can only guarantee a global minimum when the function to be minimized 
is quadratic, as L 7 (x, A) is quartic in x, the method is nonetheless efficient. It is important 
to note also that a directional minimizing algorithm employing cubic interpolation does not 
produce step sizes that would guarantee a global minimum through the conjugate gradient 
method even for a quadratic function; this feature is key in producing an AL algorithm that 



max I 0, q(x) + 
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does not easily become trapped in local minima. 

Essentially, the AL method employed seeks to minimize u = R 2 iteratively with an 
increasing penalty for overlap between spheres. This can be seen as beginning with N 
permeable spheres, squeezing them together within a spherical boundary of radius R, and 
then iteratively decreasing their permeability such that they force the boundary outward. 
The AL method also does not guarantee a global minimum of (Al), though for large enough 
values of 7 in certain problems, 7 > 70, a local saddle point of (A2) is guaranteed to 
exist [31]. For more detailed information on the augmented Lagrangian and conjugate 
gradient methods, we refer the reader to one of the many texts on the subject of nonlinear 
programming, such as [31]. 

The local minimum of (Al) found by the AL method for large enough 7 is dependent 
not only on the algorithm parameters but also strongly on the initial conditions for x and 
A. Accordingly, the procedure we use to find global minima begins with a variety of initial 
conditions for x, from disks positioned randomly according to the Poisson distribution inside 
or on the surface of a disk of a larger radius to disks positioned with centers on the sites of 
the triangular lattice. After each iteration of the AL method, a subset of the disk centers 
arranged as the local minimum in x are moved radially inward by a random amount not 
exceeding their initial distance from the origin. Immediately after, these same disks are 
rotated through a random angle such that they are no farther than distance unity from 
their previous position. This 'shuffled' configuration is used as the initial conditions for the 
first Nd variables of x in the next AL iteration, while the final variable of x is set to the 
greatest squared distance from the origin to any of the N disk centers. The initial A at 
every iteration is set to zero to help keep the AL algorithm from becoming stuck in a local 
minimum. 

In our trials, generally no more than 50 iterative AL and shuffling steps were necessary 
to find a DLP optimal packing, even for large numbers of disks (for small numbers of disks, 
as few as 1 or 10 steps was often sufficient). To support the conjecture that the minima 
found are indeed global minima, we repeated each procedure of 50 iterations as many as 20 
times, changing the governing parameters of the AL algorithm and employing the different 
initial conditions discussed above. For the vast majority of iV in the packings presented in 
this work, the same least local minima were found in all or the majority of the 20 repetitions 
performed. 
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Appendix B: Characteristics of the DLP optimal packings presented 

The following table lists the putative R min (N), calculated with accuracy to at least 1CT 6 
diameter units, for each packing presented. For comparison with the triangular lattice, 
alongside each R min (N) is listed the smallest radius triangular lattice shell enclosing the 
centers of at least N disks (not including the central disk). Point group symmetries are 
additionally listed, as determined also with accuracy to at least 1CT 6 diameter units. 
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